In kmcd39/divM: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = T,
                      fig.width=6, fig.height=5,
                      out.width = "100%")

# setup ws ----------------------------------------------------------------
library(sf)
library(tidyverse)
# rm(list=ls())
devtools::load_all()
# source(here::here("R/Generate-measures/rays/setup ray ws.R"))
load(here::here("R/Generate-measures/ray-ws.Rdata"))

# refresh crs -------------------------------------------------------------

# sometimes gets unbundled from object when moving across systems
st_crs(hwys) <- "+proj=lcc +lon_0=-90 +lat_1=33 +lat_2=45"
st_crs(plc) <- "+proj=lcc +lon_0=-90 +lat_1=33 +lat_2=45"


# get cleaned interstate hwy plan ---------------------------------------------------

plan.dir <- "/scratch/gpfs/km31/other-local-data/1947plan/most addl cleans/"
hwyplan <- list.files(plan.dir, pattern = "\\.shp", full.names = T)

#hwyplan = st_read("~/R/shapefiles/1947plan/addl cleans 2/cleaner-hwy-plan.shp")
hwyplan <- st_read(hwyplan)
# mimic structure of NHPN built hwys using plan data
# makes using the same functions easier.
hwyplan$SIGNT1 = "plan"
hwyplan$SIGN1 = paste0(hwyplan$SIGNT1,hwyplan$id)

# spatial clean hwy data --------------------------------------------------
# FIX NOLA & houston (check)
'library(mapedit)
library(leaflet)
mapview(hwyplan)
hpc <- hwyplan %>% filter(id != 161)
hpc %>% mapview()
hpc2 <- mapedit::editFeatures(hpc)
hpc2 %>% mapview()
intpl <- st_transform(hwyplan
                      ,st_crs(plc))

library(lwgeom)
sum(!st_is_valid(intpl))
intpl <- st_set_precision(intpl, 1000)
intpl <- st_make_valid(intpl)
st_write(hpc2, "~/R/shapefiles/1947plan/addl cleans 2/cleaner-hwy-plan.shp")'

# associate with place geoids ---------------------------------------------
intpl <- st_transform(hwyplan
                      ,st_crs(plc))

Illustrating hwy plan rays

There is possibility of a buffer to handle noise for low-reso hwy plan differently.

#plc <- st_buffer(plc, 8046.72)

Looking at selected places

To look at how well interstate plan seems to generate. Contrasts built- and planned- interstates in Philly at the end

plc.ids <- plc$geoid
names(plc.ids) <- plc$plc.name

sd <- Count.rays(
  plc.ids[grepl("San Diego", names(plc.ids))]
           ,intpl
           ,plc
           ,always.include = "plan"
           ,min.segment.length = 10
           ,include.map = T
)

sd$map

nola <- Count.rays(
  plc.ids[grepl("New Orleans", names(plc.ids))]
  ,intpl
  ,plc
  ,always.include = "plan"
  ,min.segment.length = 10
  ,include.map = T
)

nola$map

# philly interstates plan vs built
phl.planned <- 
  Count.rays(
  plc.ids[grepl("Philadelphia", names(plc.ids))]
  ,intpl
  ,plc
  ,always.include = "plan"
  ,min.segment.length = 10
  ,include.map = T
)

# ~built interstates with NHPN data
phl.built <- 
  Count.rays(
  plc.ids[grepl("Philadelphia", names(plc.ids))]
  ,hwys
  ,plc
  ,always.include = "I"
  ,min.segment.length = 10
  ,include.map = T
)

phl.planned$map
phl.built$map